function h2dinit(params)

  tic;

  Nx = params.Nx;
  Ny = params.Ny;
  Dx = params.Dx;
  Dy = params.Dy;
  


  % ****************************************************************
  %  compute derivative kernels
  % ****************************************************************
  % can upgrade on this -- trefethen pg 6

  % the derivative code assumes Nx = Ny so far

  global d1x d1y d2x d2y;
  global gd1x gd1y gd2x gd2y;

  d1x = zeros(Nx, Nx);
  d1y = zeros(Ny, Ny);
  d2x = zeros(Nx, Nx);
  d2y = zeros(Ny, Ny);

  % x derivatives
  for x = 1:Nx
    for y=1:Nx
      s1 = 0;
      s2 = 0;
      for k = -Nx/2 : Nx/2
        s1 = s1 + 2. * pi * i / Nx / Nx * k * exp(2 * pi * i * k * (x - y) / Nx);
        %s2 = s2 + (2. * pi * i / Nx * k)^2 / Nx * exp(2 * pi * i * k * (x - y) / Nx);
      end
      d1x(x,y) = real(s1)/Dx;
      %d2x(x,y) = real(s2)/Dx/Dx;
    end
  end

  % y derivatives
  for x = 1:Ny
    for y=1:Ny
      s1 = 0;
      s2 = 0;
      for k = -Ny/2 : Ny/2
        s1 = s1 + 2. * pi * i / Ny / Ny * k * exp(2 * pi * i * k * (x - y) / Ny);
        %s2 = s2 + (2. * pi * i / Ny * k)^2 / Ny * exp(2 * pi * i * k * (x - y) / Ny);
      end
      d1y(x,y) = real(s1)/Dy;
      %d2y(x,y) = real(s2)/Dy/Dy;
    end
  end

  global chebM1;
  global chebM0;

  N = params.Nz;
  [D, x] = cheb(N-1);

  ii  = 1:N-1;
  chebM1 = inv(D(ii,ii));

  ii  = 2:N;
  chebM0 = inv(D(ii,ii));

  % the other method
  % h = 2*pi/GRID_SIZE;
  % col1 = [0 .5*(-1).^(1:GRID_SIZE-1).*cot((1:GRID_SIZE-1)*h/2)]  * h/DX;
  % max(max(d1 - toeplitz(col1, col1([1 GRID_SIZE:-1:2]))));
  % col2 = [(-pi^2/(3*h*h)-1/6) .5*(-1).^(0:GRID_SIZE-2).*(csc((1:GRID_SIZE-1)*h/2).^2)]  * h*h/(DX*DX);
  % d2 = toeplitz(col2);

  %gd1x = gpuArray(d1x);
  %gd1y = gpuArray(d1y);

  tder = toc;
  disp(sprintf('Initializing derivative kernels, %.2f sec', toc));
end
